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Abstract 

We study the Master equation with time-dependent coefficients, a hnear kinetic 
equation for the Markov chains or for the monomolecular chemical kinetics. For the 
solution of this equation a path summation formula is proved. This formula repre- 
sents the solution as a sum of solutions for simple kinetic schemes (kinetic paths), 
which are available in explicit analytical form. The relaxation rate is studied and 
a family of estimates for the relaxation time and the ergodicity coefficient is devel- 
oped. To calculate the estimates we introduce the multi-sheeted extensions of the 
initial kinetics. This approach allows us to exploit the internal ("micro") structure 
of the extended kinetics without perturbation of the base kinetics. 
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1 Introduction 



1.1 The problem 



First-order kinetics form the simplest and well-studied class of kinetic systems. It 
includes the continuous-time Markov chains |l|2j (the Master Equation [3] ) , kinetics 
of monomolecular and pseudomonomolecular reactions [3], provides a natural lan- 
guage for description of fluxes in networks and has many other applications, from 
physics and chemistry to biology, engineering, sociology, and even political science. 
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At the same time, the first-order kinetics are very fundamental and provide the 
background for kinetic description of most of nonhnear systems: we ahnost always 
start from the Master Equation (it may be very high-dimensional) and then reduce 
the description to a lower level but with nonlinear kinetics. 

For the description of the first order kinetics we select the species-concentration lan- 
guage of chemical kinetics, which is completely equivalent to the states-probabilities 
language of the Markov chains theory and is a bit more flexible in the normalization 
choice: the sum of concentration could be any positive number, while for the Markov 
chains we have to introduce special "incomplete states" . 

The first-order kinetic system is weakly ergodic if it allows the only conservation law: 
the sum of concentration. Such a system forgets its initial condition: the distance 
between any two trajectories with the same value of the conservation law tends 
to zero when time goes to infinity. Among all possible distances, the h distance 
(11^ Ik = X^i plays a special role: it decreases monotonically in time for any 
first order kinetic system. Further in this paper, we use the li norm on the space of 
concentrations. 

Straightforward analysis of the relaxation rate for a linear system includes com- 
putation of the spectrum of the operator of the shift in time. For an autonomous 
system, we have to find the "slowest" nonzero eigenvalue of the kinetic (generator) 
matrix. For a system with time-dependent coefficients, we have to solve the linear 
differential equations for the fundamental operator (the shift in time). After that, 
we have to analyze the spectrum of this operator. Beyond the simplest particular 
cases there exist no analytical formulas for such calculations. 

Nevertheless, there exists the method for evaluation of the contraction rate for the 
first order kinetics, based on the analysis of transition graph. For this evaluation, 
we need to solve kinetic equations for some irreversible acyclic subsystems which we 
call the kinetic paths (jlOp . These kinetic paths are combined from simple fragments 
of the initial kinetic systems. For such systems, it is trivial to solve the kinetic 
equations in quadratures even if the coefficients are time-dependent. The explicit 
recurrent formulas for these solutions are given ()12p . 

We construct the explicit formula for the solution of the kinetic equation for an 
arbitrary system with time-dependent coefficients by the summation of solutions of 
an infinite number of kinetic paths psp . 

On the basis of this summation formula we produce a representation of the li con- 
traction rate for weakly ergodic systems (j23p . Because of monotonicity, any partial 
sum of this formula gives an estimate for this contraction. 

To calculate the estimates we introduce the multi-sheeted extensions of the initial 
kinetics. Such a multi-sheeted extension is a larger Markov chain together with a 
projection of its (the larger) state space on the initial state space and the following 
property: the projection of the extended random walk is a random walk for the 
initial chain (Section 14. 2p . 
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This approach allows us to exploit the internal ("micro") structure of the extended 
kinetics without perturbation of the base kinetics. 

It is difficult to find, who invented the kinetic path approach. We have used it in 
1980s [5j, but consider this idea as a scientific "folklore". 

In this paper we study the backgrounds of the kinetic path methods. This return to 
backgrounds is inspired, in particular, by the series of work |6|7] . where the kinetic 
path summation formula was introduced (independently, on another material and 
with different argumentation) and applied to analysis of large stochastic systems. 
The method was compared to the kinetic Gillespie algorithm [8j and on model 
systems it was demonstrated [7| that for ensembles of rare trajectories far from 
equilibrium, the path sampling method performs better. 

For the linear chains of reversible semi-Markovian processes with nearest neighbors 
hopping, the path summation formula was developed with counting all possible 
trajectories in Laplace space [9j. Higher order propagators and the first passage 
time were also evaluated. This problem statement was inspired, in particular, by 
the evolving field of single molecules (for more detail see |10j). 

The idea of kinetic path with selection of the dominant paths gives an effective 
generalization of the limiting step approximation in chemical kinetics |ll|12j . 



2 Basic Notions 



Let us recall the basic facts about the first-order kinetics. We consider a general 
network of linear reactions. This network is represented as a directed graph (digraph) 
( |13|14j ): vertices correspond to components Ai {i = 1,2,... ,n, edges correspond 
to reactions Ai — )• Aj {i ^ j). For the set of vertices we use notation A, and 
for the set of edges notation £. For each vertex, Ai ^ A, a, positive real variable 
Cj (concentration) is defined. Each reaction Ai — >■ Aj is represented by a pair of 
numbers («, j), i 7^ j. For each reaction Ai — t- Aj a nonnegative continuous bounded 
function, the reaction rate coefficient (the variable "rate constant") kji{t) > is 
given. To follow the standard notation of the matrix multiplication, the order of 
indexes in kji is always inverse with respect to reaction: it is /cj^j, where the arrow 
shows the direction of the reaction. The kinetic equations have the form 



or in the vector form: c = K{t)c. The quantities Cj are concentrations of Ai and c is 
a vector of concentrations. We don't assume any special relation between constants, 
and consider them as independent quantities. 

For each t, the matrix of kinetic coefficients K has the following properties: 




(1) 
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• non-diagonal elements of K are non-negative; 

• diagonal elements of K are non-positive; 

• elements in each column of K have zero sum. 



This family of matrices coincides with the family of generators of finite Markov 
chains in continuous time ( |1|2] ). 

A linear conservation law is a linear function defined on the concentrations 6(c) = 
^■ftjCj, whose value is preserved by the dynamics ([T]). Equation ([T]) always has a 
linear conservation law: b^{c) = = const. 

Another important and simple property of this equation is the preservation of pos- 
itivity for the solution of ([1]) c{t): if Cj(to) > for all i then Ci{ti) > for ti > to- 

For many technical reasons it is useful to discuss not only positive solutions to ([1]) 
and further we do not automatically assume that Cj > 0. 

The time shift operator which transforms c(to) into c{t) is U{t, to). This is a column- 
stochastic matrix: 

Uij{t,to)>0, '^Uij{t,to) = I {t>to). 

i 

This matrix satisfies the equation: 

^^l^^ = KU{tM) or '^ = Y^{k,,{t)u,i-k,^{t)uu) (2) 

i 

with initial conditions U{tQ,tQ) = 1, where 1 is the unit operator {uij{to,to) = 6ij). 

Every stochastic in column operator [/ is a contraction in the li norm on the in- 
variant hyperplanes Yli — const. It is sufficient to study the restriction of U on 
the invariant subspace {x \ J2i = 

\\Ux\\ < S\\x\\ if ^^a^j = 

i 

for some 5 < 1. The minimum of such 6 is 6u, the norm of the operator U restricted 
to its invariant subspace {x \ Yli^i — 0}- One of the definitions of weak ergodicity 
is 5 < 1 [15]. The unit ball of the li norm restricted to the subspace {x \ Xj = 0} 
is a polyhedron with vertices 

g'' = l{e'-e^), 1 + 3 . (3) 

where e* are the standard basis vectors in M": = Jjfc, ^jfc is the Kronecker delta. 
For a norm with the polyhedral unit ball, the norm of the operator U is 

max||;7(t;)|| , 
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where V is the set of vertices of the unit ball. Therefore, for a ball with vertices ([3]) 
6u = \\U\\ = l;maxy^ \uki - Ukj\ < I . (4) 

k 

This is a half of the maximum of the li distances between columns of U . The 
ergodicity coefficient, £u = 1 — djj , is zero for a matrix with unit norm 6u = ^ and 
one if U transforms any two vectors with the same sum of coordinates in one vector 
{6u = 0). 

The contraction coefficient (jj]) is a norm of operator and therefore has a "sub- 
multiplicative" property: for two stochastic in column operators U, W the coefficient 
6uw could be estimated through a product of the coefficients 

^uw < SuSw ■ (5) 

We will systematically use this property in such a way. In many estimates we find 
an upper border 1 > 5{t) > 5(7(ti+T,ti)7 ^2 > ^i- In such a case, 5(7(ti+T,ti) ^ 
exponentially with r ^ oo. Nevertheless, the estimate 5{t) may originally have 
a positive limit 6{t) doo > when r — >■ oo. In this situation we can use 5{t) 
for bounded r < ti and for r > ti exploit the multiplicative estimate ([5]). The 
moment ri may be defined, for example, by maximization of the negative Lyapunov 
exponent: 

Ti = arg max < > . (6) 

T>0 [ T ) 

For a system with external fluxes Ili{t) the kinetic equation has the form 



^ = Y.^k^Jit)cj - kjiit)ci) + Hiit) . (7) 
j 

The Duhamel integral gives for this system with initial condition c(to): 
(t) = U{t,to)c{to) + fu{t,T)U{T) dr , 



ft 

'to 

where n(r) is the vector of fluxes with components nj(r). 



In particular, for stochastic in column operators U{t, to) this formula gives: an iden- 
tity for the linear conservation law 

ft 



Y,Ci{t) = Y,Ci{to)+ [ Y.U,{T)dT , (8) 

and an inequality for the l\ norm 

<||[/(t,to)c(to)||+ r5]||n(r)||dr<||c(to)||+ T ^ ||n(r)|| dr . (9) 

Jtn- J tn 



*0 j *0 . 

We need the last formula for the estimation of contraction coefficients when the 
vector c(t) is not positive. 
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3 Kinetic Paths 



Two vertices are called adjacent if they share a common edge. A directed path is a 
sequence of adjacent edges where each step goes in direction of an edge. A vertex 
A is reachable from a vertex B, if there exists a directed path from B to A. 

Formally, a path in a reaction graph is any finite sequence of indexes (a multiindex) 
/ = {ii, ^2, • • • iq} ((? > 1, 1 < ij < n) such that (ifc, ife+i) € £" for all A; = 1, . . . , g — 1 
(i.e. there exists a reaction Ai^^ — )■ Aj^,^^). The number of the vertices |/| in the path 
/ may be any natural number (including 1), and any vertex Ai can be included in 
the path / several times. If g = 1 then we call the one-vertex path I degenerated. 
There is a natural order on the set of paths: J >/ if J is continuation of /, i.e. 
/ = {zi, i2, • • • iq} and J = {ii,i2, ■ ■ ■ iq, ■ ■ ■}■ In this order, the antecedent element 
(or the parent) for each / is I~, the path which we produce from / by deletion of 
the last step. With this definition of parents /~, the set of the paths with a given 
start point is a rooted tree. 

Definition 1 For each path I = {ii,i2, ■ ■ ■ iq} we define an auxiliary set of reaction, 
the kinetic path Pj: 

2(J2) q{lg) ^^QS^ 



The vertices Bj^.^^^ of the kinetic path (jlOp are auxiliary components. Each of them 
is determined by the path multiindex / and the position in the path /. There is a 
correspondence between the auxiliary component Bj^.^s^ and the component A^ of 
the original network. The coefficient is a sum of the reaction rate coefficients for 
all outgoing reactions from the vertex Ai of the original network, and the coefficient 
K-7 is this sum without the term which corresponds to the reaction Ai Aj-. 

f^i — ^ ^ ^/i) '^jj = ^ ^ kii . 
I, Ij^i I, 



A quantity, the concentration b^^-^y corresponds to any vertex of the kinetic path 
BI^-^^ and a kinetic equation of the standard form can be written for this path. The 
end vertex, Bj^^- y plays a special role in the further consideration and we use the 
special notations: ij = iq, Aj = Ai^, qj = b^g(^i^y i^i is the reaction rate coefficient of 
the last outgoing reactions in (jlOp (the last vertical arrow) and kj is the reaction 
rate coefficient of the last incoming reaction in (jlOp (the last horizontal arrow). 

We use for the incoming flux for the terminal vertex of the kinetic path (jlOp 
and Pf for the outgoing flux for this vertex. 

Let us consider the set Zi of all paths with the same start point ii and the solutions 
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of all the correspondent kinetic equations with initial conditions: 



= 1' ^m) = for / > 1 . 

For the concentrations of the terminal vertices this self-consistent set of initial con- 
ditions gives the infinite chain (or, to be more precise, the tree) of simple kinetic 
equations for the set of variables I £ Ii: 

?i = -Ki(t)<;-i, ^/ = -Ki{t)qi + ki{t)c;i- , (11) 

where index 1 corresponds to the degenerated path which consists of one vertex 
(the start point only) and corresponds to Ai^ . 

This simple chain of equations with initial conditions, <?i(to) = 1 and ^/(to) = for 
\I\ > 1, has a recurrent representation of solution: 



?i(i) = exp ^- j^^ Ki(r) dr^ , 

?/(i) = ^ exp (^-^ K/(r)dr^ kiie)^i^{e) dO 



(12) 



The analogues of the Kirchhoff rules from the theory of electric or hydraulic circuits 
are useful for outgoing flux of a path J G Xi and for incoming fluxes of the paths 
which / are the one-step continuations of this path (i.e. I~ = J): 

i-=J 

Let us rewrite this formula as a relation between the outgoing flux PJ from the last 
vertex of J and incoming fluxes for the last vertices of paths / (/^ = J): 

pj= E (14) 

/, i-=J 



The Kirchhoff rule (jl4p together with the kinetic equation for given initial conditions 
immediately implies the following summation formula. 

Theorem 1 Let us consider the solution to the initial kinetic equations ([7]) with 
the initial conditions Cj(to) = <5jii- Then 

/ell, ii=j 



Proof. To prove this formula let us prove that the sum from the right hand side 
(i) exists (ii) satisfies the initial kinetic equations ([T]) and (iii) satisfies the selected 
initial conditions. 
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Convergence of the series with positive terms follows from the boundedness of the 
set of the partial sums, which follows from the Kirchhoff rules. According to them, 

because Ti consists of the paths with the selected initial point ii only. 
The sum 

/gll, ii=j 

satisfies the kinetic equation ([T]). Indeed, let Iij = {I € Ii \ ii = j} be the set of all 
paths from ii to j. Let us find the set of all paths of the form {/~ | / G Iij}- This 
set (we call it Z^) consists of all paths to all points which are connected to Aj by 
a reaction: 

X,- = l) lu. 

From this identity and the chain of the kinetic equations (jlip we get immediately 
that 

(16) 

The coincidence of the initial conditions for q and Ci is obvious. Hence, because of 
the uniqueness theorem for equations ([T]) we proved that Ci = Ci. □ 

It is convenient to reformulate Theorem [T] in the terms of the fundamental operator 
U {t, to)- The ith column of U{t, Iq) is a solution of ([1]) Cj{t) = Uji{t, to) {j = 1, . . . ,n) 
with initial conditions Cj{to) = 6ij. Therefore, we have proved the following theorem. 
Let Iij be the set of all paths with the initial vertex Ai and the end vertex Aj and 
<^/(t) be the solutions of the chain pT]) for ii = i with initial conditions: ?i(to) = 1 
and ?7(to) = for |/| > 1. 

Theorem 2 

Uji{t,to) = Y^i{t) . □ (17) 



Remark 1. It is important that all the terms in the sum (jlTh are non-negative, 
and any partial sum gives the approximation to Uji{t,tQ) from below. 

Remark 2. If the kinetic coefficients are constant then the Laplace transform gives 
a very simple representation for solution to the chain (jlip (see also computations 
in [9|6j ). The kinetic path / (jlOp is a sequence of elementary links 

k' k 



The transfer function Wi^{p) for the link (fTSj) is the ratio of the output Laplace 
Transform to the input Laplace Transform for the equation. Let the input be a 
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function Xj,, (t) and the output be Yi^{t) = > where bi^{t) is the solution to 
equation 

kl = -K-hK + 1 K = -l^irK + kirir-l^iri't) > 1) 

with zero initial conditions. The Laplace transform gives 

for a link (jlSp and for the whole path (jlOp we get 

Wi = —^\\^^^. (19) 

(compare, for example, to formula (9) in [6j). It is worth to mention commutativity 
of this product: it does not change after a permutation of internal links. For the 
infinite chain pip with initial conditions, ?i(0) = 1 and ?/(0) = for |/| > 1, the 
Laplace transformation of solutions is 

= Wi (20) 



4 Evaluation of Ergodicity Coefficient 



4-1 Preliminaries: Weak Ergodicity and Annihilation Formula 



4.1.1 Geometric Criterion of Weak Ergodicity 

In this Subsection, let us consider a reaction kinetic system ([T]) with constant coef- 
ficients kji > for G £■ 

A set E is positively invariant with respect to the kinetic equations ([T|), if any 
solution c{t) that starts in E at time to (c(to) S E) belongs to S for t > to (c(t) S E 
if t > to)- It is straightforward to check that the standard simplex S = {c | Cj > 
0, J2i c« = 1} is a positively invariant set for kinetic equation ([T]): just check that if 
Cj = for some i, and all Cj > then q > 0. This simple fact immediately implies 
the following properties of K: 

• All eigenvalues A of ii' have non-positive real parts, Re\ < 0, because solutions 
cannot leave S in positive time; 

• If ReX = then A = 0, because the intersection of S with any plane is a polygon, 
and a polygon cannot be invariant with respect to rotations to sufficiently small 
angles; 

• The Jordan cell of K that corresponds to the zero eigenvalue is diagonal - because 
all solutions should be bounded in E for positive time. 
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• The shift in time operator exp(i^t) is a contraction in the li norm for t > 0: there 
exists such a monotonically decreasing (non-increasing) function 5{t) (t > 0, 
< 6{t) < 1, that for any two solutions of ^ c{t),c'{t) G S 

Y,h{t)-C^^{t)\<S{t)Y,\cM-cm\■ (21) 
i i 

Moreover, if for c{t),c'{t) € S the values of all linear conservation laws coincide then 
\ci{t) — c'^{t) \ — )■ monotonically when t — )• oo. 

The first-order kinetic system is weakly ergodic if it allows only the conservation 
law: the sum of concentration. Such a system forgets its initial condition: distance 
between any two trajectories with the same value of the conservation law tends to 
zero when time goes to infinity. 

The difference between weakly ergodic and ergodic systems is in obligatory existence 
of a strictly positive stationary distribution: for an ergodic system, in addition, a 
strictly positive steady state exists: Kc = and all Cj > 0. Examples of weakly 
ergodic but not ergodic systems: a chain of reactions Ai —?■ A2 ■ ■ ■ An and 
symmetric random walk on an infinite lattice. 

The weak ergodicity of the network follows from its topological properties. 

Theorem 3 The following properties are equivalent ( and each one of them can he 
used as an alternative definition of weak ergodicity): 

(1) There exists a unique independent linear conservation law for kinetic equations 
(this is lP{c) = Cj = const J. 

(2) For any normalized initial state c(0) (b^{c) = 1) there exists a limit state 

c* = lim exp(/s:t) c(0) 

that is the same for all normalized initial conditions: For all c, 

lim exp(Kt)c = 6°(c)c*. 

(3) For each two vertices Ai, Aj {i ^ j) we can find such a vertex Aj. that is 
reachable both from Ai and from Aj . This means that the following structure 
exists: 

Ai^ ...^ Ak^ ...^ Aj . (22) 
One of the paths can he degenerated: it may he i = k or j = k. 

(4) For t > operator exp{Kt) is a strong contraction in the invariant suhspace 
^iCi = in the h norm: || exp(i('t)a;|| < 6{t) < 1, the function 6{t) > is 
strictly monotonic and 6{t) when t — t- 00 □. 
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The proof of this theorem could be extracted from detailed books about Markov 
chains and networks ( |1I17| ). In its present form it was published in [5j with explicit 
estimations of the ergodicity coefficients. 

Let us demonstrate how to prove the geometric criterion of weak ergodicity, the 
equivalence 1 3. 

Let us assume that there are several linearly independent conservation laws, linear 
functionals b^{c), b^{c), . . . ,b^{c),m > 1. The linear transform c i— >■ (6^(c), . . . , 6"^(c)) 
maps the standard simplex S„ in M" (cj > 0, q = 1) onto a polyhedron D C M™. 
Because of linear independence of the system b^{c), 6^(c), . . . , 6'"(c), m > 1, this D 
has nonempty interior. Hence, it has no less than m + 1 vertices wi, . . . ,Wq, q > m. 

The preimage of every point x € -D in S„ is a positively invariant subset with 
respect to kinetic equations because the standard simplex is positively invariant 
and the functionals 6*(c) are the conservation laws. In particular, preimage of every 
vertex Wq is a positively invariant face of Fq; FqD Fr = if 7^ r. 

Each vertex Vi of the standard simplex corresponds to a component Ai: at this 
vertex q = 1 and other Cj = there. Let the vertices from Fq correspond to the 
components which form a set Sq; Sq H Sj- = ^ q ^ r. 

For any Ai G Sg and every reaction Ai — )• Aj the component Aj also belongs to Sq 
because Fq is positively invariant and a solution to kinetic equations cannot leave 
this face. Therefore, if g 7^ r, G Sq and Aj E Sr then there is no such vertex 
that is reachable both from Ai and from Aj. We proved the implication 3 =^ 1. 

Now, let us assume that the statement 3 is wrong and there exist two such com- 
ponents Ai and Aj that no components are reachable both from Ai and Aj. Let Si 
and Sj be the sets of components reachable from Ai and Aj (including themselves), 
respectively; Si fl Sj = 0. 

For every concentration vector c G M" a limit exists c*(c) = limj_j.oo exp(i^t) c 
(because all eigenvalues of K have non-positive real part and the Jordan cell of 
K that corresponds to the zero eigenvalue is diagonal). The operator c i-^- c*(c) is 
linear operator in M". Let us define two linear conservation laws: 

b\c)= C:(C), V{C)= ^ c;(c) . 

These functionals are linearly independent because for a vector c with coordinates 
Cr = Sri we get 6*(c) = 1, 5^ (c) = and for a vector c with coordinates Cr = 5rj we 
get 6*(c) = 0, V{c) = 1. Hence, the system has at least two linearly independent 
linear conservation laws. Therefore, 1^3. 

4.1.2 Annihilation Formula 

Let us return to general time-dependent kinetic equations ([I]). 
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In this Section, we find an exact expression for the contraction coefficients 6{t, Iq) for 
the time evohition operator U{t, to) in li norm on the invariant subspace {x \ Xi = 
0}. The unit /i-ball in this subspace is a polyhedron with vertices g^^ = ^(e' — e^), 
where e, are the standard basic vectors in M" The contraction coefficient of an 
operator U is its norm on that subspace ([H), this is half of the maximum of the li 
distances between columns of U. 

The kinetic path summation formula ()17p estimates the matrix elements of C/(t, to) 
from below, but this does not give the possibility to evaluate the difference between 
these elements. To use the summation formula efficiently, we need another expression 
for the contraction coefficient. 

The ith column of U{t, to) is a solution of the kinetic equations ([1]) Cj{t) = Uji{t, to) 
(j = 1, . . . ,n) with initial conditions Cj{to) = 5ij. For each j let us introduce the 
incoming flux for the vertex Aj in this solution: 

mt)=Y,k,g{t)cg{t) 

(the upper index indicates the number of column in U{t,to), the lower index corre- 
sponds to the number of vertex Aj). 

Formula ([H) for the contraction coefficient gives 

5{t,to) = lmax\\U{t,to){e' - e^)\\ . 

U{t,to){e' — e^) is a solution to the kinetic equation ([T]) with initial conditions 
Ci{to) = 1, Cj{to) = —1 and Cq{to) = for q ^ This is the difference between 
two solutions, c^(t) = Ugi{t,to) and c~{t) = Uqj{t,to). Let us use the notation 

G'^{t) = lu{t,to){e'-e^) . 

For each q we define 

n^= E kgi{c+-c^), n-= J2 kqi{q-c+), n±>o. 

l,c+>c- '.c+<c- 

The decrease in the li norm of c^{t) — c~{t) can be represented as an annihilation 
of a flux n^(t) with an equal amount of concentration c~^{t) —c~{t) from the vertex 
Aq by the following rules: 

(1) If Cq = c^{t) — c^{t) > then the flux 11" annihilates with an equal amount 
of positive concentration stored at vertex Aq (Fig. [Ta|) ; 

(2) If Cq = Cq{t) — Cq{t) < then the flux annihilates with an equal amount 
of negative concentration stored at vertex Aq (Fig. [Tbll ; 



12 




c" 




(a) c > 0, the nega- 
tive flux annihilates 



(b) c < 0, the posi- 
tive flux annihilates 

Fig. 1. Annihilation of fluxes. 



(c) c = 0, the mini- 
mal flux annihilates 



(3) If Cq = c^{t) — (t) = then the flux min{n+, 11^ } annihilates with the equal 
amount from the opposite flux (Fig. [Tc]) . 

Let us summarize these rules in one formula: 
Proposition 1 



dt' 



q, C^>Cq 



q, c+<c. 



E min{n+(t),n-(t)} . □ 



(23) 



q, C^=Cq 



Immediately from (j23p we obtain the following integral formula 



1 = / I E n-(r)+ E n+(r) 

\q, C^>Cq q, C^<Cq 

+ E min{n+(r),n-(r)}) I dr 



(24) 



q, c ' =c„ 



The annihilation formula gives us a better understanding of the nature of contraction 
but is not fully constructive because it uses fluxes from solutions to the initial kinetic 
equation ([T]). 



4.2 Multi-Sheeted Extensions of Kinetic System 



Let us introduce a multi-sheeted extension of a kinetic system. 

Definition 2 The vertices 0/ a multi-sheeted extension of the system ([1]) are AxK 

where K is a finite or countable set. An individual vertex is {Ai,l) (I € K). The 
corresponding concentration is C(^i^i)- The reaction rate constant for {Ai,l) — )■ {Aj,r) 
is k(j,r){i,i) ^ 0- This system is a multi-sheeted extension of the initial system if an 
identity holds: 

E hj,r)(i,i) = kji for all I . (25) 
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Fig. 2. Redirection of a reaction from one sheet to another with preservation of the 
base kinetics. The redirected reaction is highhghted by bold. 

This means that the flux from each vertex is distributed between sheets, but the 
sum through sheets is the same as for the initial system. We call the kinetic behavior 
of the sum q = Yli the base kinetics. 

A simple proposition is important for further consideration. 

Proposition 2 If C(^i^i^ (t) is a solution to the extended multi-sheeted system then 

C^{t) = Y,Ci^,l)(t) (26) 
I 

is a solution to the initial system and 

Elc(M)WI>El^^(*)l • (27) 

il i 

(Here we do not assume positivity of all ci). □ 

Formula (j25l) allows us to redirect reactions from one sheet to another (Fig. [2]) 
without any change of the base kinetics. In the next section we show how to use 
this possibility for effective calculations. 

Formula ()26p means that kinetics of the extended system in projection on the initial 
space is the base kinetics: the components {Ai,l) are projected in Ai the projected 
vector of concentrations is q = C(jj) and the projected kinetics is given by the 
initial Master equation with the projected coefficients kji = ^{j,r){i,i)- "Recharg- 
ing" is any change of the non- negative extended coefficients k(^j^r){i,i) which does not 
change the projected coefficients. 

The key role in the further estimates plays formula (j27p . We will apply this formula 
to the solutions with the zero sums of coordinates, they are differences between the 
normalized positive solutions. 

4-3 Fluxes and Mixers 

In this Subsection, we present the system of estimates for the contraction coefficient. 
The main idea is based on the following property which can be used as an alternative 
definition of weak ergodicity (Theorem [3]) : For each two vertices Ai, Aj (i ^ j) we 
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can find a vertex Ag that is reachable both from Ai and from Aj. This means that 
the following structure exists: 

y . . . y ^ ' ' ' ^ * 

One of the paths can be degenerated: it may he i = q or j = q. The positive flux 
from Ai meets the negative flux from Aj at point Ag and one of them annihilates 
with the equal amount of the concentration of opposite sign. 

Let us generalize this construction. Let us fix three different vertices: Ai (the "pos- 
itive source"), Aj (the "negative source") and Ag (the "mixing point"). The degen- 
erated case q = i OT q = j we discuss separately. Let be such a system of vertices 
that Ai G S~^, Ag ^ and there exists an oriented path in U {Ag} from Ai 
to Ag. Analogously, let S~ be such a system of vertices that Aj £ S~ , Ag ^ S~ 
and there exists an oriented path in S~ U {Ag} from Aj to Ag. We assume that 

s+nS' = 0. 

With each subset of vertices S we associate a kinetic system (subsystem): for A^ (z S 

n 

Cr = ^ ^ k^lCl — ^ ^ kprCr . (28) 
I, Ai£S, r^l p=l 

In this subsystem, we retain all the outgoing reaction for Ar G S and delete the 
reactions which lead to vertices in S from "abroad" . 

The flux Ug from to Ag is 

= ^ ^ kgrCrit) , 

r, AreS+ 

where Cr{t) is a component of the solution of (128p for S = S*^ with initial conditions 
Cr(io) = ^ri- Analogously, we define the flux 

Lltj = ^ ^ kgrCr(t) , 
r, Ar&S- 

where Cr{t) is a component of the solution of (I28p for S = S~ with initial conditions 
Cr(*o) = (^rj- Decrease of the norm ||G*-'(t)|| is estimated by the following theorem. 

The system S~^,S~,Ag we call a mixer, that is a device for mixing. An elementary 
mixer consists of two kinetic paths Ai — )• ... —^Ag ^ ... ^ Aj ()22p with the 
corespondent outgoing reactions: 

y . . . > Ai^ i 



where ii = i, ir = q, ir+i = j- 



A. 



(29) 
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Fig. 3. A mixer: two subsystems, (includes Ai) and S~ (includes Aj). There 
may be outgoing reactions from but all incoming reactions to from outside 
are deleted. A mixing point Aq and two fluxes, positive from (marked by dark 
color) and negative from S~ , meet at the mixing point. 

The degenerated elementary mixer consists of one kinetic path: 

(30) 

K. — K. — K • 

n'2 »2»3 '■r 

where ii = i, ir = j- 



Theorem 4 



G'^{t)\\ < 1 - / min{n+,ns} dt . (31) 



Proof. To prove this theorem let us organize a 4-sheeted extension of the initial 
kinetic system as it is demonstrated in Fig. [3l Subsystems including the positive 
source (initial concentration +1) and the negative source (initial concentration — 1) 
belong to level 0. Reactions from to Aq are redirected to the sheet /, reactions 
from S'^ to other vertices, which do not belong to S^, go to sheet +1, reactions 
from to other vertices, which do not belong to , go to sheet —1. The incoming 
flux to the sheet / is 11^ — 11^ . 

Let us introduce the following notations: 

n 

ApeS+ q=l 
n 
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r=l 

We consider solution to the kinetic equations for the multi-sheeted system with 
initial conditions: C(j o)(io) = 1; C(j_o)(^o) = ~1 and all other concentrations are 
equal to zero at time to- In this case, some of the signs of concentrations are known 
for t >tQ due to the organization of the redirection of reactions (Fig. [3]) : 

C(p^o) ^ G , C(p 0) < for Ap £ S~ , 

0) = for Ap^S+US- , (32) 

C(q,l) > 0, % _i) < . 

Let us use dH) for with the sheet +1 and for S~ with the sheet —1. We get 
immediately 

dCt „^ dC, 



— ^ = n+ — ^ = n~ (33) 



Analogously, we can use Q for the sheet / and get 

<|n+-n-|. (34) 



dCf 



For the norm of the base vector of concentration c the inequality holds (Proposi- 
tion [2]): 

||c|| <Cj + Cs+Cf. 



Finally, we combine this inequality with (I33p . (I34p and get 

||c(t)|| < 2 -2 / min{n+(r),n_5(T)} dr □ 

Jtn 



For the degenerate case the path from Ai goes directly to Aj (or inverse), let us 
assume that there is a subsystem , Ai € , the mixing point Aq coincides with 
Aj and the flux is 

r, ArG5+ 

where Cr{t) is a component of the solution of ()28p for 5 = with initial conditions 

Cr{to) = Sri- 

Theorem 5 

/•min{t,ii } 

||G^^-(t)||<l- / n+(T)dT, (35) 

J to 

where Kj = kpj and ti is a solution to equation 

t 

nj(r) exp(— Kj(t — r)) dr = exp(— Kjt) . (36) 



h 
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Proof. This theorem is also proved by the construction of the appropriate multi- 
sheeted extension of the kinetic system. For the degenerated case we need only two 
additional sheets: subsystem including the positive source Ai (initial concentra- 
tion +1) and the negative source Aj (initial concentration —1) belong to level 0. 
Reactions from to other vertices, which do not coincide with Aj, go to sheet +1, 
reactions from Aj to other vertices go to sheet —1. The concentration of ^(j,o) is 

Hj,o)i't) = [ (r) ex.p{-Kj{t - r)) dr - exp{-Kjt) . 
Jto 

Let us introduce the following notation: 

n 

= ^(P.O) + m ^(9.1) ' 

ApeS+ q=l 

n 
9=1 



For t <ti concentrations C(j o)(*) ^^(9,-1) negative, hence 

and for the norm of the correspondent solution for the base system we get the 
inequality 

/•min{t,ti } 

||c(t)|| < 2-2 / n+(r) dr □ (38) 

Jto 



The kinetic path summation formula gives us a family of estimates of Hg from 
below. For each pair i,j we can find the best of available estimates of ||G*-'(t)|| (the 
smallest one for various choices of Aq and subsets S^) and then among all pairs 
of i,j we should choose the "most pessimistic" evaluation of ||G*-^ (i)|| (the biggest 
one). It will give the evaluation of the contraction coefficient from above. 



5 Simple example: Irreversible Cycle 

Let us demonstrate all results for a simple kinetic system, a simple irreversible cycle: 

Ai^A2^... ^ An ^ Ai (39) 

All ki > and are constant in time. For enumeration of A^ we use the standard 
cyclic order (modn): A^^j = Aj. 
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The kinetic equations for this system are: c = Kc or 





Cl 




d 


C2 

















~ki ... kn 
ki -k2 ... 

-kn 
The characteristic equation for this system is 





Cl 




C2 




_ c„ _ 



(40) 



i=l 



i=l 



One eigenvalue for matrix K is, obviously, A = 0, the correspondent left eigenvector 
is the linear conservation law li = (1, 1, ... , 1). The right eigenvector for this A is 
the steady state ri = j_ • • • ' at)'^ (normalized for /iri = 1). Other n — 1 

roots of the characteristic equations have strictly negative real parts, ReXi < 
{i > 1) and, in general, cannot be found explicitly. For a given eigenvalue A, the 
eigenvectors have a simple structure: 



hi 



+1 



Ar 



X + kj 



-7—, V'Ai-l 
Ki 



X + kj 
ki 



(41) 



With the normalization condition: for eigenvalues A, A': ^a^a' = ^xx'^ that is 1 for 
A = A' and for A A'. 

Two limit cases allow explicit analysis of eigenvalues and eigenvectors of K: 



(1) Systems with limiting steps: one constant is much smaller than others, let it 
be kn, kn<^ki, {i = l,...,n- 1); 



(2) Fully symmetric systems, ki = k2 



kn- 



For systems with limiting steps {kn <C ki, {i = 1, . . . ,?7- — 1)) the eigenvalues are 
close to —ki, . . . , —kn-i and the relaxation is limited by the second constant, the 
next to the minimal one (detailed analysis is provided in |11|12] ). 



For a symmetric system {ki = k2 — ... — n,„ — n,;, ui±c cig,ci±vaiuco ci±c. xvg 
/cexp ^^^^ — 1 for q = 1, . . . ,n. There are n distinct eigenvalues, one of them. 



kn 



k), the eigenvalues are: Xq 



k 



cos 



/ 2niq \ 



\ n I 



1 



. Let us 



An = 0, the other have negative real part: ReXq 

further take k = 1 for this system (include k into dimensionless time). For the 
left and right eigenvectors (|41|) we have two waves moving in opposite directions, 

Iqj+i = ^gi^xpf^^j, Vqj-i = rqjexpl^^Y We can take with respect to the 
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normalization condition, IqVp = Sqpi 



l>exp ^^^^ ,exp ^2^^^ , . . . ,exp ^(n - 1)^^^ 

1 / / 2TTiq\ f 27riq\ 

- 1, exp , exp -2 , . . . , exp 

n \ \ n I V n I 



27r«g\ 

-(n- 1) 

n / 



(42) 



For constant coefficients, the operator of sliift in time from to ti depends only on 
t = ti — to: U{ti,to) = U{t) = expKt. We can use (|^ and write 



C/(t) = ^exp(A,t)|r,)(/,|, 



g=l 



{U{t))js = ^exp(Agt)rgj^q3 



(43) 



9=1 



1 " 

- Vexp 



9=1 



t I cos 1 

n 



cos (s — J j h r sm 

' n n 



This explicit formula allows us to compute all the necessary quantities including 



the contraction coefficient (5, 



Uit) 



Now, let us produce the approximate formula for the same symmetric system by 
mixers. First of all, let us represent the solution for the cycle by the path summation 
formula. With the convention of cyclic enumeration, the set of paths Xj started at 
Ai is the sequence 



1i= { 



hi 



Ai — y Ai+i > , 



Ai % ^,+1 ^ A,, 



^1+2 



. . . 



> . 



(44) 



This sequence of paths corresponds to the multi-sheeted representation presented 
in Fig. m First, we consider a infinite series of the copies of the cycle. Each vertex of 
the extended system is numerated by two indexes: {Ai,l), i = 1, 2, . . . , n (modn), I = 
1,2,3,... is a natural number. The reaction rate constants for copies are the same 
as for the initial systems: k(^j^j.^(^ij^ = kjiSri- This extended system obviously satisfies 
the definition of the multi-sheeted extension of the cycle and in its projection on 
the base we always have the kinetics of the cycle. 

Let us select one number i E {1, . . . ,n} and recharge the reactions: we annulate 
the "horizontal" reaction rate constant for {Ai,l) — >■ (j4j+i,/), ^(j+i, = 0, and 
instead of this reaction take the reaction between levels, {Ai,l) — > (^j+i,/ + 1): 
= ki+ii (see Fig. H]). This is also a multi-sheeted extension of the 
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r 




Fig. 4. Multi-sheeted representation of the path summation formula for a cycle (|46p : 
a cycle (the base) is represented by an semi-infinite helix produced by redirecting 
of reactions between sheets. 



cycle. Formula (126p for this multi-sheeted system allows us to use integration of the 
infinite acyclic system (represented by the spiral in Fig. [H)) instead of integration 
of the finite cyclic base system. 

Now, let us put all ki = 1. For systems with constant coefficients we use initial time 
moment Iq = 0. For the set of paths li started at Ai the solution to the chain (jlip 
with the initial conditions ft (to) = 1 and = for |/| > 1 is 

^/(i) = (^^--' • (45) 
Obviously, X^/gj. = 1- For concentration of Ag, formula (fT7|) gives 

where dij is the length of the shortest oriented path from Ai to Aj (here the length 
is the number of reactions and the trivial path from Ai to Ai has the length zero). 

For every two vertices Ai, Aj we have only two mixers and both are degenerated: 

Ai A- ^i+i ^ . . . Aj A-, length j — i mod n and Aj A- vlj+i A- . . . A- 
length i — j mod n. 

Let us select one mixer Ai ^ A2 ■ ■ ■ Aj A- for analysis. Initial conditions are: 
ci = 1, Cj = — 1 and other concentrations are equal to zero. 

For this auxiliary chain with given initial conditions 



{p-l)\ 



e-* (p=l,...,j-l). 



-e-* I 1 



-1 X (47) 



{j - 1) 
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The estimate ([35]) ||G*-^ (t)|| < 1 — /q n^(r) dr is valid until cj changes its sign. 
Hence, for t we have a boundary < (j — 1)!. The Stirling formula gives a 
convenient estimate: 

i-i 



t^-' < v/2vr(i - 1) (^)' < (i - 1)! 

7-1 1 

t < ti = (27r(j - l))2(j-i) . 

e 



(48) 



Even a simpler estimate is t < {j — l)/e. If t satisfies one of these inequalities then 
concentration Cj is negative and we can use the estimate (j35p . 

For this example, 

n+(t) = c,_i(t) = T^e-* , rn+(r) dr = 1 - e-* X; ^ , 

mm{dji,d,j}-l [2J 

p=0 ^' p=0 ^' 

where [^] is the integer part of n/2. For t > this estimate gives ||G*''(t)|| < 1 
and 5(7(t) < 1 because X]p=o ff ^ '^^'^ estimate (|19|) on an interval 

[0, ti], for example, on [0, ^^] - Intersection of these intervals for all i,j, i ^ j is [0, j] 
(j > 2). On this interval, the estimate (j49p is valid for all For extension of such 
an estimate for i > ^ the submultiplicative property ([5]) can be used. 



6 Ergodicity Boundary and Limitation of Ergodicity 



In this Section we consider a reaction kinetic system ([T]) with constant coefficients 
kji > for (z, j) e £. 

Let us sort the values of kinetic parameters in decreasing order: k(i-^ > /c(2) > . . . > 
A;(„). The number in parenthesis is the number of value in this order. Each of the 
constants is a reaction rate constant kij for some i,j (and may be for several 
of them if values of these constants coincide). Let us also suppose that the network 
is weakly ergodic. We say that /c^,,), 1 < r < n is the ergodicity boundary [18] if 
the network of reactions with parameters ki,k2, ■ ■ ■ ,kr is weakly ergodic, but the 
network with parameters ki,k2, ■ ■ ■ , kr-i is not. In other words, when eliminating 
reactions in decreasing order of their characteristic times, starting with the slowest 
one, the ergodicity boundary is the constant of the first reaction whose elimination 
breaks the ergodicity of the reaction digraph. 

Let A4ij {i 7^ j) be a set of elementary mixers (f29|) . (pO]) between given Ai, Aj. For 
each M € A4ij we can find a cutting reaction rate constant, cut a/: 

cutM = uim{ki2ii ki^i^_^ , , . . . , } for ([29]) ; ^^^^ 

cutM = mm{ki2ij^,. . . , ki^i^_^] for ([30]) . 
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Let us eliminate reactions in increasing orders of their constants (i.e. in decreasing 
order of their characteristic times), starting with the smallest one. To cut all ele- 
mentary mixers between Ai,Aj {i ^ j), it is necessary and sufficient to eliminate 
all kpg < cutM for all M £ Aiij. Therefore, for every pair Aj, Aj {i / j) we can also 
introduce a cutting constant: 

cutji = max cutM • 

To destroy the weak ergodicity of the network J\f we have to cut al least one pair 
Ai,Aj (i ^ j). The result can be formulates as the following theorem. 

Theorem 6 Theorem 6 The ergodicity boundary of a network J\f is the following 
constant: 

cut AT = min cutjj . □ 

This boundary is a minimum (in pairs Ai,Aj) of maxima (in mixers M G Mij) of 
minima (in constants). 

Kinetic equations for elementary mixers ()29p . (j30p allow explicit analytic solutions. 
Nevertheless, explicit estimates in terms of cutting constants can be also useful. 

Let for an elementary mixer M ()29p km be the maximal sum of constants of outgoing 
reactions: 

KM = raax{Kip \ p = ii,i2,. . . ,ir+i}, = ^ kps , 

or for a degenerated elementary mixer M ()30p 

KM = max{Kip \ p = ii,i2,. ■ ■ ,ir} ■ 



Let us substitute all the constant for horizontal arrows in the elementary mixer M 
(f29l) . (pOl) hy k = cutAf , and all the constants for vertical arrows {i ^ ir) hy k — k, 
where k = km- This change decreases the fluxes H^. 

To find the estimate we have to solve the kinetic equation for a simple uniform 
kinetic path: 



Al — ^ A2 — ^ ... — ^ As 



K—k 



K—k 



K—k 



(51) 



Similar to the simple cycle ()47p . we find 

{kt)P~^ 



" (p-1)! 
the only difference is in exponents. 



exp(-Kt) (p= l,...,s) , 



(52) 
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For the elementary mixers (j29p . ()30p this formula gives 



and the estimates from Theorems I4I5I (|31 p . (j35p become simple analytical expressions 
after substitution of by their estimates from below. 

Let us find an universal estimate from below for ti. It is 



k + K 



Indeed, in the degenerated elementary mixer ()30p on the way from Ai to Aj there 
exists at least one reaction with reaction rate constant k: A^ .... The integral 
flux through this reaction during the time interval [0, t] is 



/ kcrir) dr > / U+{t) dr . 
Jo Jo 



The last inequality holds because all the flux in the mixer should go through the 
reaction Ar ^ . . . before it enters the last vertex. On the other hand, kcr{T) dr < 
Jp* k exp(— A;r) dr (the last integral corresponds to the case when all the concentra- 
tion is collected at the initial moment at Ar and goes only through the reaction 
Ar — 7> . . .). Therefore, 



/ n+(r) dr < 1 - exp(-/cr) . 
Jo 



From the condition (I36p we find the estimate for ti from below: ti > ri, where ri is 
solution to 

1 — exp(— A:r) = exp(— Kr) . 

We use convexity of exponential functions and substitute them in this equation by 
linear approximation at point r = 0: exp(— x) > 1 — x {x > 0); this gives us the 
estimate of ri from below: ri < = -j^^- 

For t G [0,??], kt<l and 

{ktl [ktl [kty 
0! 1! ^ ■■■ ^ r\ 

For each mixer M we introduce the length of mixer (Im = max{r — 2,1 — 1} for 
(|29p and cIm = r — 2 for ()30p . In these notations, each mixer M € A4ij gives the 
estimate: for t € [0, 1?^/] 

WG^'m < 1 - f cntM^^^^^^^^eM-i^MT) dr , (53) 

where 



cuIm + Km 
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For each pair i,j [i ^ j) we can select the "critical" elementary mixer M G A4ij with 
cutAf = cutjj and put dij = cIm, i^ij = km- If there are several critical elementary 
mixers then we select one with minimal cIm , if there are several such a mixers with 
minimal dM then we select one with minimal km- In this notation we have 



for t G [0, "i^ij], where 



< 1 







* (cut.jr)'^'^^ 
cutjj — — exp(— KjjTj cir 



{di,y- 



(54) 



1 



Finally, for the whole network 

cut^r = min {cutjj}, dj\f = max{dij}, Kj\f = max{Kjj}, "dj^ = - 

i,3,i^3 * J:«7^i «J:«7^i CUt^T + Hj\f 

and for the contraction coefficient 5{t) ([2T]) we obtain the estimate 



5{t) <1 



* (cutATT)'^-^ . . , 

— 73-"Ti — exp(-K_ArT) dr 







(%)! 



1 ~ — □ — exp(-KArt) 



p=0 



(55) 



for t € [0, For t outside this interval, the submultiplicative property ([S]) should 
be used. 



7 Discussion 



The kinetic path summation formula together with the multi-sheeted extension of 
kinetics provide us with a factory of estimates. It is difficult to find, who invented 
this approach. 

The analysis of kinetic paths with selection of the most important (dominant) paths 
allowed us to extract dominant systems from kinetic equations |ll|12j . A robust pro- 
cedure for simplification of biochemical networks was created [19]. This approach 
was developed into unified framework for hybrid simplifications of Markov models of 
multiscale stochastic gene networks dynamics [20] . Dominant subsystems were ana- 
lyzed for dynamical models of microRNA action on the protein translation process 

The multi-sheeted extension of kinetics provides us with a simple and useful tech- 
nique for estimation of relaxation processes in Master equation. This method intro- 
duces an internal "microstructure" in the first order kinetic systems. The kinetic 
path summation formula is a particular case of the formula (126p (Proposition [2]) . 
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Indeed, let us construct the following multi-sheeted extension of the Master equa- 
tion. The set of components is A x IC, where K, = {0} U fCi and )Ci is the set of 
all kinetic paths / with lengths |/| > 1 (non-degenerated paths). The connections 
between sheets (redirected reactions) are: 

Ai^_j- ^ Aijj instead of ^ Ai^j~ . 

According to this rule, the reaction that continues the path /~ to the path / is 
redirected and goes from the sheet /~ to the sheet I. For a degenerated /~, we take 
A^^_ j- = Ai^_^Q, this means that all paths start on the zero sheet, and all reactions 
from this sheet lead to other sheets: Ai — )• Aj transforms into Ai^o — >■ Aj^^^ jy, where 
{i,j} is a path of the length 2. Formula (|26p for this multi-sheeted structure coincide 
with the kinetic path summation formula (jl7p (Theorem [2]) for initial conditions 
Cifi = 1 and other cq- /) = 0. 

This multi-sheeted extension may be considered as a generalization of the Bethe 
lattices introduced by H. Bethe in 1935 [22]. For example, if in the initial graph of 
reactions each vertex has the same number of outgoing edges then the constructed 
multi-sheeted extension can be considered as a bundle of the Bethe lattices, each 
of them starts from one point of the zeroth sheet. For each starting point, ^(i,o) 
the corresponding Bethe lattice represents the "Green function" Uji{t,to) for given 
i and for all possible j. 

We produced the kinetic path summation formula for time-dependent kinetic equa- 
tions and applied this formula for evaluation of the ergodicity coefficient. The eval- 
uation of the contraction coefficient in the li norm is the main tool for studying of 
the relaxation in time-dependent Markov processes since the seminal works of R. 
Dobrushin |15j . 

Another important context of this study is the analysis of the eigenvalues of the 
stochastic matrices |23|24j and, especially the analysis of these eigenvalues for ma- 
trices with specified graph |25|26j . In chemical kinetics, evaluation of the eigenvalues 
through kinetic constants was given in series of work by V.Cheresiz and G. Yablon- 
skii [27128] . 

Various estimates of eigenvalues of K could be produced from the estimates of 
contraction ()3ip . (135^ . The simplest one follows from (I55p : 

«e(A) < < . (56, 

v 

Several problems should be resolved to make the use of the path summation formula 
more effective. Perhaps, the most important of them was mentioned in the comment 
[29]). The amount of the kinetic path needed for accurate estimate of the solution 
grows quickly in time for a sufficiently complex system. Hence, we need either special 
tricks for the analysis of path sampling or special asymptotic formulas for long paths 
instead of exact solutions. 

Another possible approach to this problem is in the use of more complex exactly 
solvable systems instead of paths. The set of reactions is solvable, if there exists 
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a linear transformation of coordinates c i— )■ a such that kinetic equations in new 
coordinates for all values of reaction constants have the triangle form: 

^ = fi{ai,a2, ..-ai). (57) 
at 

The algorithm for the analysis of reaction network solvability was developed in [5] 
(see also [TT] ) . The simplest examples of solvable networks give acyclic graphs (reac- 
tion trees) and pairs of mutually inverse reactions. It may be possible to decompose 
the complex system of transitions into a sequence of solvable systems. 
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